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Strain relaxations of a p(l x 2) ordered oxygen layer on W(llO) are measured as a function of 
temperature across the disordering transition using low-energy electron diffraction. The measured 
strains approach values of 0.027 in the [110] and —0.053 in the [001] direction. On the basis of the 
measured strain relaxations, we give quantitative information on temperature-dependent surface 
stress using the results of ab initio calculations. From the surface formation energy for different 
strains, determined by first-principles calculations, we estimate that surface stress changes from 
— 1.1 for the ordered phase to —0.2 N/m for the disordered one along [110], and from 5.1 to 3.4 N/m 
along [001]. Moreover, our observation that the strains scale inversely with domain size confirms 
that the strain relaxation takes place at the domain boundaries. 



I. INTRODUCTION 

Surface stress is a quantity important for the under- 
standing of surface processes like reconstruction, inter- 
facial mixing, segregation, critical film thickness and 
self-organization at solid surfacesii"— Regarding the lat- 
ter, recent studies have demonstrated that the stress- 
induced patterns occurring at the mesoscopic scale on 
metal and semiconductor surfaces^"— show a strong de- 
pendence on temperature,— suggesting that the change of 
surface stress with temperature can play a major role in 
determining their properties. In spite of the wide inter- 
est in these and other stress-driven phenomena, there are 
very few studies concerning the temperature dependence 
of surface stress 

On the theoretical side, density-functional theory 
(DFT) is probably the most accurate method to study 
surface elastic properties; however its results apply to 
Kiii Non-zero temperatures can be investigated us- 
ing model potentials and atomistic simulations, but in 
the literature there are few examples applied to surface 
elasticity as a function of temperatures^^ More stud- 
ies can be found on a related quantity, the surface free 
energy^i^^— which, however, does not give information on 
surface stresses without knowledge of its dependence on 
the surface strain. On the experimental side, the crystal- 
bending method was used to measure stress changes due 
to adsorbed species at different temperatures Indi- 
rectly connected to surface stress, there are some mea- 
surements of temperature-dependent strain using x-ray 
diffraction. However, a truly temperature-dependent 
study of a given system is still lacking. 

Here, we aim at estimating surface stress change as a 
function of temperature across the order-disorder tran- 
sition of the p{l X 2)-O/W(110) structure, using an al- 
ternative approach which combines low-energy electron 
diffraction (LEED) and ab initio calculations. In addi- 



tion to thermal expansion, the presence of the disordering 
is crucial in determining the surface stress changes, sim- 
ilar to the phenomenon of surface meltingJ^ This stress 
change is also of particular importance for the recently 
reported high-temperature stress- induced pattern forma- 
tion of Pd on W— , due to the profound effect of O ad- 
sorption on the anisotropy of Pd stripes 

Oxygen on W(llO) is a well-studied model system 
characterized by a series of ordered phases as a func- 
tion of oxygen coverage J^^i^ At 0.5 ML coverage O is 
ordered in a (1 x 2) structure, which consists of doubly- 
spaced close-packed rows of O in the (111) directions. 
The (1 X 2) phase contains eight equivalent domains gen- 
erated by translations and rotation (by 109.5°) of the O 
lattice on the triply coordinated adsorption sitesJ^ The 
p(l X 2) O domains disorder at about 700 K. Bucholz and 
Lagally^S describe the disordering by shrinking of the or- 
dered domains with increasing temperature. At the same 
time, the surface stress is expected to be relieved at the 
domain boundariesJ^ 

In this paper, we use LEED to detect the temperature- 
induced lattice changes in the O adlayer on W, and DFT 
calculations for their quantitative analysis in terms of 
surface stress. More precisely, we estimate the surface 
stress by using surface formation energy as a function of 
strain, calculated from first principles. We evaluate the 
surface stress at the values of strain for which the sur- 
face layer is found to be maximally relaxed (well above 
the disordering temperature). We focus on the behavior 
in a wide temperature range, and we show that the av- 
erage change in atomic spacing of the surface layer for 
increasing temperature is a direct consequence of stress 
relaxation due to disorder on the surface. The analy- 
sis of lattice changes measured by LEED is done in the 
same way as in Rcf. [2ll . Kinematic low energy electron 
diffraction, although similar to x-ray diffraction in the in- 
formation sought, has the advantage of being inherently 
surface sensitive, even more so in the case of O/W(110) 
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FIG. 1: (a) LEEM image of the partially (1x2) oxygen-covered 
W(llO). Image diameter is 2 microns. The curved lines cor- 
respond to W atomic steps, (b) LEED pattern acquired from 
a 2 micron region at 35 eV electron energy, (c) and (d) Dark 
field LEEM images of the same region acquired using the half- 
order spots marked in panel (b). The bright/dark contrast 
corresponds to the 1x2 and 2x1 domains. 



with the presence of half-order diffraction spots. 

The paper is organized as follows: in section II we 
present the measured strain relaxations as a function of 
temperature, obtained by LEED. Section III summarizes 
our ab initio calculations of surface formation energy as 
a function of strain. Combining the two, we extract es- 
timates of surface stress at high temperature, which is 
discussed in section IV. Conclusions are given in section 
V. 



II. EXPERIMENTAL RESULTS 

The LEED measurements were made with the spec- 
troscopic photo-emission and low energy electron micro- 
scope SPELEEM at Elettra (Italy) .-^ The instrument is a 
hybrid electron/x-ray microscope, which allows chemical 
and structural imaging combined with micro-spot LEED 
(^-LEED)i^i2^ In the latter mode, a LEED pattern was 
obtained from a region of 2 fim diameter. The trans- 
fer width of /i-LEED was determined to be 10 nm from 
measurements on a step-free region of the clean W(llO) 
surface. 

Cleaning of the W(llO) single crystal was carried out 
by annealing at 1400 K in 2 x 10~^ mbar oxygen, followed 
by high temperature flashes in ultrahigh vacuum to re- 
move oxygen. The surface cleanliness was confirmed by 



the absence of tungsten carbide islands and a sharp (1x1) 
LEED pattern. Using low-energy electron microscopy 
(LEEM) , regions with micron-sized terraces were chosen 
in order to exclude broadening of diffraction spots due to 
the presence of step bunches. 

The (lx2)-ordered 0.5 ML oxygen-covered surface 
was prepared by exposing the clean W(llO) surface to 
molecular oxygen at 450 K. 0.5 ML coverage was as- 
signed to the maximum intensity of the (1x2) diffrac- 
tion spots, which was obtained at an exposure of about 
4.6 Langmuirsj^ After O exposure, we performed a short 
annealing at 1400 K under ultra-high vacuum conditions, 
which resulted in large domains with a sharp LEED pat- 
tern. 

Figure [T^ shows a LEEM image of the (1x2)- 
O/W(110) surface acquired below the disordering of oxy- 
gen, at around 400 K. The corresponding /z-LEED pat- 
tern is displayed in Fig. [T|d. The very large (1x2) and 
(2x1) domains are visible in the dark field images shown 
in Fig. [1]: and[T]ii^ The domain sizes reach micron scale 
with boundaries consisting of narrow regions with oppo- 
site domain orientation between two rotational domains 
showing no obvious directional anisotropy. 

It should be noted that the presence or absence of 
translational domain boundaries within a region of sin- 
gle orientation (1x2 or 2x1) cannot be concluded from 
the images displayed in Fig. [TJ However, the high inten- 
sity of the half-order diffraction spots confirms the pres- 
ence of extended single domains. Nevertheless, fine lines 
and granularity in a single rotational domain in the high 
resolution images could point to the presence of defects 
(translational domain boundaries, oxygen vacancies, etc) 
within an otherwise perfect oxygen lattice. 

This surface is used to monitor the strain relaxations 
at high temperature. The disordering of the (1x2) or- 
der as a function of temperature can be followed via the 
intensity and width of the half-order diffraction spots dis- 
played in Fig. [5^. The decrease of intensity and the in- 
crease of spot widths are documented in the early litera- 
ture^SiS^ and can be understood as large domains break- 
ing up into smaller ones with increasing temperature. 
This is also reflected in the shape of the diffraction spot 
profiles, which transforms from Gaussian to Lorentzian 
at about 600 K. Importantly, diffuse half-order spots can 
still be discerned up to above 1100 K. This means that 
the O atoms preserve the local (1x2) order (in small 
clusters) even though the long-range order is lost at high 
temperature. 

In addition to the intensity and width, we can also 
measure the half-order peak positions as a function of 
temperature. The reciprocal lattice distances allow us 
to evaluate the average distance of O atoms within the 
(1 X 2)-ordered regions. The data analysis is carried out in 
the same way as in Ref. [111. The resulting strains, defined 
with respect to the underlying W lattice, in the two crys- 
tallographic directions are displayed in Fig. [2}d. Clearly, 
the (1x2) O unit cell, which respects the W periodicity 
at low temperatures, on average expands along [110] and 
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FIG. 2: (a) Peak intensity (solid line) and width of the half- 
order diffraction spots, measured at electron energy of 28 eV. 
(b) Measured strain, derived from the reciprocal distances of 
the half-order peak positions, as a function of temperature 
along two directions. Solid lines are the least-squares fit, in- 
dicated to guide the eye. 



shrinks along [001] as the temperature increases and dis- 
order sets in. At the highest temperatures the measured 
strains tend to reach values of about eroiaxod = 0.027 in 
[110] and eroiaxed = -0.053 in [001] directions at 1100 K. 

The measured values of the strain relaxation include a 
contribution from the thermal expansion of the W crys- 
tal. However, within the temperature range of interest, 
the change in lattice spacing due to this effect is less than 
OA %^ which is comparable to the smallest measured 
surface strains. 

Interestingly, the strains in both directions scale as in- 
verse power laws of the average domain size as shown in 
Fig. [3] Domain sizes were extracted from the FWHM of 
half-order diffraction spots as displayed in Fig. [T]i tak- 
ing proper account of the instrumental broadening. The 
fit functions in the log-log plots in Fig. [3] are described 
by £ cx L^, and correspond to powers of p = —1.09 and 
— 1.45 for [110] and [001] directions, respectively. This 
observation is in line with the strain relaxations taking 



place at the domain boundaries. Moreover, the changes 
in strain as a function of O coverage at fixed temperature 
(450 K) are also included in Fig.[3]for comparison)^ The 
important differences will be discussed in the following 
sections. 

Note that the two main sources of error in the ex- 
periments arc the limited instrument transfer width and 
the changes in the angular alignment function of 

temperature. The former makes the determination of 
domain sizes very difficult when the domains are larger 
than the transfer width (hence the broadening of the 
diffraction spots is dominated by the instrument func- 
tion). However, considering that the transfer width is 
about 10 nm, this should be a negligible effect within the 
range displayed in Fig. [31 On the other hand, variations 
in angular alignment may introduce a small uncertainty 
in the spot positions, which might translate into slight 
deviations from the power law behavior. In order to re- 
duce the possibility of such an artifact, we have averaged 
over equivalent lattice vectors in obtaining the reported 
strains. 



III. AB INITIO CALCULATIONS 

In this study, wc use ab initio calculations to determine 
the surface formation energy, 7, for different strains. The 
strain derivative is used to evaluate the surface stress at 
the measured maximally-relaxed strain ereiaxed of the O 
layer. 

To calculate 7 we perform DFT pseudopotential cal- 
culations in a plane-wave basis by using the PWscf 
codci^ The local-density approximation (LDA) in the 
Perdew-Zunger parametrizationS^ is used for the ex- 
change and correlation functional. We utilize a sym- 
metric slab with 11 layers to simulate the 0/W system, 
with 9 vacuum layers. We employ Vanderbilt ultra-soft 
pseudopotentials^^ generated using the 2s^2p'* atomic 
configuration of O and the 5s^5p^5d*6s^ configuration 
of W. The core-cutoff radii for O were: rg^p = 1.6, 
rd — 1.4 a.u. and for W: rs^p = 2.2, = 2.4 a.u. We 
use the equilibrium W bulk lattice constant of 3.14 A. 
Our plane-wave-basis kinetic-energy cutoff is 35 Ry for 
the wave functions and 350 Ry for the charge density. 
We use a 26 X 13 X 1 Monkhorst-Pack^il mesh for the 
k-points sampling in the Brillouin zone of the (1 x 2) 
surface. 

7 is evaluated from the expression: 
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where A is the strain-dependent surface area, -Bsiab 
stands for the total energy of the whole slab, N^y for 
the number of W atoms, for the W bulk energy per 
atom. No for the number of O atoms and E02— is the 
energy of the O2 molecule. Esiab is evaluated for several 
different homogeneous uniaxial [110] and [001] strains on 
the slab and E^^ is recalculated for each strain. Wc note 
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FIG. 3: Changes of strain along (a) [110] and (b) [001] as a function of the average dimension of the (1 x 2) domains. The 
squares and circles are the data extracted from the LEED measurements, while the solid lines are fits to the power-law functions. 
Temperatures corresponding to 20, 50 and 100 A domain sizes are noted above. Strains obtained in Ref. [2l| as a function of 
oxygen coverage at fixed temperature are included as dashed lines for comparison. 



that the constant energy E02 does not contribute to the 
derivative of energy and, therefore, is not relevant to the 
discussion of surface stress. 

On the basis of the calculated 7 for various strains 
and a fitted polynomial curve, the surface stress, is 
obtained by using)^ 



1 djAj) 



1,2 



(2) 



where Ci is the uniaxial strain. Figure |4] presents 7 mul- 
tiplied by the area^i for a few strains, as well as poly- 
nomial fits to the calculated points for both directions. 
We assume that the two directions can be considered in- 
dependently. In the [110] direction, a quadratic function 
fits the data very well, while for the [001] direction, due 
to higher-order terms, a better fit is obtained by a cubic 
polynomial. 

Calculated values of r are presented in Table HI In the 
first column are our p{l x 2)O/W(110) surface stress val- 
ues calculated from the curves in Fig. |4]at e = 0. The 
second column contains the same quantities, for compari- 
son, from Ref. [2l| (calculated with the expression derived 
by Nielsen and Martin2^ for the stress tensor a based 
on the Hellmann-Feynman theorem) for both directions. 
Finally, in the last column are our estimates for the av- 
erage surface stress in the high-temperature phase, with 
no long-range order {T^isO{ix2)), as obtained from the 
derivative of the energy curve in Fig. 2] at eroiaxod- 



IV. DISCUSSION 

The calculated values in TableUfor the stresses at e = 
using the two approaches, from the 7 curve and by a 
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TABLE I: Surface stresses at e = calculated from the 7 curve 
(first column) and by a direct DFT surface stress calculatioi*^ 
(second column). Last column contains surface stresses at 
Erciaxcd, estimated from the 7 curve. The values are given in 
N/m. The positive sign corresponds to the tensile stress. 



direct DFT calculation of the surface stress, give similar 
results in both directions. The differences are within the 
estimated error, 0.35 N/m,— of the calculated surface 
stress and the propagated error of the fitting parameters 
of ~ 0.20 N/m. The different slab thicknesses in the two 
calculations (15 layers for the direct stress calculation and 
11 layers for the total energy calculation of the strained 
slabs) also may contribute to the difference. The surface 
stress values at eroiaxcd (high-temperature phase) indicate 
a significant reduction with respect to the values at e = 
(0 K). 

To better understand the physical significance of the 
values listed in the last column in Table U we need to de- 
scribe the underlying assumptions. Most importantly, we 
assume that the disordering takes place predominantly 
within the O layer and that the temperature changes in 
stress are mostly due to this disorder. The basis for the 
first approximation is that W atoms do not participate 
in the (1 x 2) order to begin withi^"— Therefore, they 
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FIG. 4: Ab initio calculations of the surface formation en- 
ergy multiplied by the surface unit cell area as a function of 
strain in a) [110] and b) [001] directions. Black dashed curves 
show quadratic fit, while the solid red line is for a cubic fit. 
Full circles indicate calculated points. The vertical line cor- 
responds to the measured strain, treiaxed, of the O layer at 
high temperature. 



are expected to be disturbed little by the disordering of 
tiie O lattice, if at all. The second assumption is based 
on the small effect of the bulk thermal expansion of W. 

As we noted in the introduction, with increasing tem- 
perature, the O lattice loses long-range order by breaking 
of large domains into smaller ones. It is expected that 
the loss of order should result in a reduction of stress 
within the O layer with the stress relieved at the domain 
boundaries In the limiting case of complete disorder 
(a dense gas of O atoms) one would expect negligible 
stress in the O layer, and the remaining surface stress 
would be mostly dominated by the underlying W layer. 
On the other hand, the presence of the half-order diffrac- 
tion spots at the highest temperatures measured, how- 
ever weak and diffuse, suggests that the O layer has still 
short-range order. Therefore, the high temperature sur- 
face stress can be decomposed into contributions from 
the unmodified W layers and from the partially relaxed 



oxygen surface. 

In order to represent this situation, in our theoreti- 
cal analysis, ideally, we would like to calculate stresses 
for the strains within the oxygen layer equal to ereiax 
and keeping the bulk unstrained in the same slab. How- 
ever, such inhomogencous strain is not affordable in our 
supercell calculations. Therefore, we use slabs in which 
both the bulk and the surface part are strained to ^relaxed 
(with subsequent subtraction of the bulk contributions). 
To see how this theoretical treatment corresponds to the 
experiment in which we measure only the surface-layer 
strain with W underneath at its equilibrium lattice con- 
stant (strain zero), we examine the layer-resolved stresses 
as a function of strain. From Eq. [T] and [2l the uniaxial 
surface stress can be written as: 

^i(ei) = 5 7-[CTi,slab(ej) - iVlvCTj^wl^i)] 

= ^^[fi,surf(ej) + fi.subsurf (eO + {N^r - l)(Tj^^(ei) 

(3) 

where ai^sUh{^i) = {l/V)dEs\ah{t{)/ dti and cr?w('^0 = 
{I /V)dE^{ei) / dei are the volume-averaged stresses, V 
is the volume of the unit cell and the indices "surf" and 
"subsurf" stand for the surface O and the subsurface W 
layers, respectively. Assuming that the subsurface layer's 
proximity to the surface changes its stress-strain relation 
only by addition of a constant stress (which corresponds 
to the stress at zero strain, see the Appendix), with re- 
spect to the bulk, we are left with: 

1 V 

= --^ [^i^surf (ej) + 0-i,subsurl(e = 0)] . (4) 

Under the above assumption, our calculated stress at 
Eroiaxod has two Contributions: from the surface O layer 
at Croiaxed and from the W subsurface layers at the strain 
zeroi^ We note that our assumption for Eq. |3] is analo- 
gous to saying that the elastic constants of the subsurface 
and bulk layers are identical. A heuristic interpretation 
of Eq. [3] and U can be found in the Appendix. 

A confirmation that the clastic response of the subsur- 
face W is close to that of bulk W is obtained by inspect- 
ing the relaxed positions in our calculated 0/W slab. We 
find that the average interlayer spacing of the two top- 
most W layers below O is very close to its bulk value 
(A(i23 = —1.1 %), while the clean W surface relaxes by 
A(ii2 — —3.6 This indicates that the W surface 

when covered with O becomes much more bulk-like. We 
note that our results for A(i23 for 0/W and Ac?i2 for W 
agree weh with other theoretical results, Ac?23 = —1.3 % 
(Ref. 113) for the former and Adi2 = -3.6 % (Ref. ^ 
for the latter case. Experimentally, Ac?i2 = —3.1 % for 
W.ii 

In spite of the assumptions mentioned above, there are 
important conclusions to draw from the results. Interest- 
ingly, our results in Table U show that the surface stress 
and its anisotropy do not vanish even several hundred 
degrees above the disordering transition at about 600 K. 
This is mainly due to short-range correlations within the 
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O layer persisting at high temperatures. The other con- 
tribution, which may be equally significant, is the non- 
zero stress within the subsurface layers which do not par- 
ticipate in the disordering. A discrimination of the two 
effects needs a layer-resolved stress calculation, which is 
beyond the scope of our study. 

Another point is the confirmation that stress is re- 
lieved at the domain boundaries, which is derived from 
the power-law behavior in Fig. [3] Although this is qual- 
itatively similar to the relaxations at the boundaries of 
small p{l X 2)0 islands on W(110)r^ there are differences 
to be emphasized. The strain relaxation at the boundary 
between two distinct phases is driven by the difference 
in surface stress, as shown in Ref. [2ll However, in the 
case of the fully 0-covered surface, the mechanism of the 
stress relaxation is somewhat different. 

To verify that, we constructed heavy and light domain- 
wall configurations and compared their energy on the ba- 
sis of broken and newly-formed bonds, using the lateral 
interaction parameters for up to the third neighbors)^ 
This simple analysis yields a vacancy line (light-domain 
wall between the (1x2) and (2x1) domains) as the most 
favorable domain-wall configuration for the main crys- 
tallographic directions on W(llO). Although the energy 
difference between this and some heavier domain walls is 
smaller than 10 meV/A, this estimate does not include 
the energy obtained by strain relaxations, which is the 
largest for this type of the boundary, as it allows the 
largest atomic relaxations. Predominant vacancy- type 
domain walls would, indeed, confirm the view that the 
stress is relaxed at the domain boundaries in this sys- 
tem. 

The differences between strain relaxation at an is- 
land boundary and a domain boundary are underlined in 
Fig.[3l Most importantly, i) along [110] strain relaxation 
on fully oxygen-covered surface as a function of domain 
size with increasing temperature shows a considerably re- 
duced magnitude (but identical slope) compared to that 
of relaxation for small oxygen islands as a function of is- 
land size with increasing coverage at fixed temperature 
(dashed line in Fig. [5^), n) along [001] both magnitude 
and slope of the strain curves are drastically different 
between two data sets. From the first point we con- 
clude that along [110] the force driving the strain relax- 
ation at the boundary is reduced for the fully 0-covered 
surface, however the relaxation mechanism remains the 
same. The island boundary responds to the stress dif- 
ference between 0-covered and clean tungsten, which is 
-4.7 N/m<^ This is, indeed, considerably larger than the 
stress on fully 0-covered surface along [110] (see Table |l| 
and it points to the idea that at the domain boundaries 
the relaxation is driven by the stress itself and not a dif- 
ference of stresses. The comparison along [001] gives an 
even stronger support to this discussion. In the case of 
an O island, the boundary almost feels no force along 
[001] and therefore the relaxation is a higher-order ef- 
fect decaying rapidly with increasing island size (dashed 
line in Fig. [5}d)- However, in the fully 0-covered sur- 



face the strain scales with a power of domain size much 
closer to a 1/L power law, which is similar to the be- 
havior along [110]. This is in line with the idea that the 
domain boundary is relaxing under the influence of the 
large tensile stress along [001]. 

We emphasize that an anisotropic thermal expansion 
of the oxygen layer cannot replace surface stress as the 
driving force of the observed strain relaxations. This is 
clearly shown in the size of the oxygen unit cell as a 
function of temperature, which can be derived from the 
data in Fig. O The oxygen unit mesh expands slightly 
from 400 K to 600 K, above which it contracts almost 
linearly with increasing temperature. The contraction is 
explained by the predominantly tensile stress within the 
oxygen layer. 

Having shown qualitatively that the strain relaxations 
are due to the stress itself and not the stress difference, 
we note that also the sign of the measured strain relax- 
ation in the two directions supports this conclusion, as 
well as the ratio of the strain relaxations (large response 
in the direction of larger surface stress). However, it is 
interesting to observe that the strain relaxation along 
[110] is only two times smaller than along [001], and not 
five times, as could be expected from the stress ratio in 
the two directions. This is explained by different poten- 
tial profiles in the two directions, according to which O 
atoms move much easier along [110], than [001]4^ 

V. CONCLUSION 

In this work, we have measured the strain relaxation 
as a function of temperature across the disordering tran- 
sition of p(l X 2)O/W(110). The breaking up of large 
domains into smaller ones as a function of temperature 
can be observed by the changes of the intensity and width 
of the half-order diffraction spots. Following the changes 
of the spot distances we extracted the strain relaxation as 
a function of temperature. At T~ 1100 K the measured 
strains reach values of 0.027 along [110] and —0.053 along 
[001]. We found that the measured stresses are roughly 
proportional to the inverse domain size. By means of 
ah initio calculations, we determined the surface forma- 
tion energy as a function of strain, from which we esti- 
mated the surface stress values at the measured relaxed 
strains. 



VI. APPENDIX 

In order to give a heuristic interpretation of Eq. [3] and 
|4] in terms of energy, in Fig. [5] we show the slab energy 
decomposed in layers^ as a function of strain. Horizon- 
tal and vertical axes represent strain and energy, respec- 
tively. The layer-energy dependence is shown shifted in 
energy, for ease of viewing. The top parabola corresponds 
to the O surface, while those below stand for the subsur- 
face and bulk W layers from top to bottom. For simplic- 
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FIG. 5: A schematic of energy decomposition in layers before 
(left) and after the subtraction of the bulk energy (right). Red 
curve represents O surface, and the full black line stands for 
the W layers. The dashed line in the second layer from the 
top represents a case with slightly changed curvature of the 
W parabola. Dots indicate continuation of bulk layers, while 
fiat lines describe zero energy contribution. The horizontal 
axis represents strain and the vertical axis stands for energy 
(contributions from different layers are drawn with energy 
shift). 



be modified with respect to bulk W. This decomposition 
is the energy analog of Eq. [S] For the subsurface W we 
show two curves, one identical to the bulk curves, and 
the other one with 25 % smaller curvature. 

When we subtract, as is done to determine surface for- 
mation energy, the bulk W energy from all the W layers, 
we are left with only two non-zero contributions, from the 
surface and subsurface layers, as shown on the right side 
of the figure. In the case when the subsurface parabola 
is identical to the bulk one, the remaining contribution 
after subtraction to the surface energy is a line, which 
will give a constant stress for all strains, as written in 
Eq. m For a subsurface parabola slightly different from 
the bulk one, the resulting contribution to the surface 
formation energy has a very small curvature, yielding a 
slowly- varying stress as a function of strain. 

Finally, adding the energy of the subsurface layer to 
the surface-layer parabola, we obtain a curve compara- 
ble to Fig. m which corresponds to the surface stress of 
the O layer and also contains contribution from the W 
subsurface layer at strain zero. 
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